Analysis
Construct a network of differentially expressed genes
# Assign nodes (treatments) and edges (5% FDR DEGs)
plotMat <- de.res %>% filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25)
allTreat <- unique(plotMat$Treatment)
target <- data.frame(Treatment = allTreat,
target = c("Autophagyi", "IAPi", "HDACi", "mTORi", "ITKi", "TLR8a", "MDM2i",
"IMID", "JAKi", "XPOi", "BCL2i"))
opt$sub <- TRUE
if(opt$sub == TRUE) {
plotMat <- de.res %>% filter(p.adj < 0.1, !Treatment %in% c("Pomalidomide", "Dacinostat")
) %>%
filter(logFC > 0.25 | logFC < -0.25)
} else {
plotMat <- de.res %>% filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25)
}
allTreat <- unique(plotMat$Treatment)
nodes <- data.frame(Treatment = allTreat)
edges <- plotMat %>% select(Treatment, Symbol) %>%
mutate(Treatment = as.character(Treatment))
# Create network object
set.seed(100)
net <- network(edges, directed = TRUE)
netTab <- ggnet2(net)$data |>
mutate(Type = ifelse(label %in% allTreat, "Treatment", "Symbol"),
Treatment = ifelse(Type == "Treatment", label, NA)) |>
select(-alpha,-color,-shape,-size) |>
mutate(x = as.numeric(x), y = as.numeric(y))
# Define edges based on direction and padj of DEGs
from <- edges %>%
left_join(netTab, by = c("Treatment" = "label")) %>% dplyr::select(Treatment, x, y) %>%
dplyr::rename(x_start = x, y_start = y) |> distinct()
to <- edges %>%
left_join(netTab, by = c("Symbol" = "label")) %>% dplyr::select(Symbol,x,y) %>%
dplyr::rename(x_end = x, y_end = y) |> distinct()
edges <- edges %>% left_join(from) %>% left_join(to) %>%
left_join(plotMat %>% select(Symbol,Treatment,logFC,p.adj)) %>%
distinct() %>%
mutate(Direction = ifelse(logFC < 0,"Down","Up"),
alpha = -log10(p.adj)
#alpha = rescale(-log10(p.adj#squish(p.adj, c(1.812407e-209,1))))
)
# Create auxiliary table for annotation
n_DEG <- plotMat %>% #filter(!Treatment %in% c("Pomalidomide", "BafilomycinA1", "Dacinostat")) %>%
group_by(Treatment) %>%
filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25) %>%
dplyr::count() %>% ungroup()
netTabAnno <- netTab |>
left_join(target) |>
left_join(n_DEG) |>
mutate(Treatment = factor(Treatment, levels = allTreat))
p <- ggplot() +
geom_curve(data = edges,
aes(x = x_start, xend = x_end, y = y_start, yend = y_end, col = Direction, alpha = p.adj),
linewidth = 0.1, curvature = 0.05) +
geom_point(data = netTabAnno[netTabAnno$Type == "Symbol",],
aes(x = x, y = y),
size = 0.5, shape = 21, color = "grey60", fill = "grey95", stroke = 0.75, alpha = 0.6) +
geom_point(data = netTabAnno[netTabAnno$Type == "Treatment",],
aes(x = x, y = y, fill = Treatment, size = n), shape = 21, stroke = 0.75) +
# scale_fill_manual(values = col_TX) +
scale_size_continuous(name = "Number of DEGs",
breaks = seq(0, max(n_DEG$n), by = 1000), limits = c(0, max(n_DEG$n)), range = c(0.5,7.5)) +
scale_color_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
scale_alpha_continuous() +
geom_label_repel(data = netTabAnno |> filter(Type == "Treatment"),
aes(x = x, y = y, label = Treatment, fill = Treatment, point.size = rescale(n, to = c(0.5,7.5))),
max.overlaps = Inf, vjust = "bottom", nudge_y = -0.001, fontface = "bold", alpha = 0.7, #direction = "y",
show.legend = FALSE, size = 5) +
guides(color = guide_legend(ncol = 2, override.aes = list(linewidth=2.5)),
size = guide_legend(ncol = 2),
alpha = guide_legend(override.aes = list(linewidth=2.5)),
fill = guide_legend(override.aes = list(size=3.5))) +
xlab("") + ylab("") + ggtitle("Network of Differentially Expressed Genes") +
lgd +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
panel.border = element_rect(linewidth = 1, colour = NA, fill = "transparent")) +
scale_fill_manual(values = c("Venetoclax" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000",
"Ibrutinib" = "#283593", "Everolimus" = "#43A047", "Selinexor" = "#FFD45F",
"BafilomycinA1" = "#F8BBD0", "Motolimod" = "#AB47BC", "Dacinostat" = "#F44336",
"Dacinostat" = "#64B5F6", "Pomalidomide" = "#80DEEA", "Ruxolitinib" = "#4DB6AC"))
p
if(opt$sub == TRUE) {
ggsave(plot = p, filename = paste0(opt$plot, "deg_subnetwork.png"),
height = 20, width = 25, units = "cm")
} else {
ggsave(plot = p, filename = paste0(opt$plot, "deg_network.png"),
height = 20, width = 25, units = "cm")
}
Volcano plots
pList <- lapply(c("Motolimod"), function(x) {
if(x %in% c("Birinapant", "Motolimod", "Everolimus", "Ibrutinib", "BaflimoycinA1")) {
geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
} else if(x %in% c("Nutlin-3a", "Selinexor")) {
geneSelec <- path.info %>% filter(pathway %in% c("TP53up")) %>% select(gene) %>% unlist() %>% as.vector()
geneSelec <- c("TP53", geneSelec)
} else if(x %in% c("Everolimus")) {
geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
}
volTab <- de.res %>% filter(Treatment == x)
## Define limits
minLim <- -1 * max(abs(de.res$logFC))
maxLim <- max(abs(de.res$logFC))
thresh <- 0.25
volTab$diffCol <- NA
volTab[volTab$logFC > 0 & volTab$p.adj < 0.1 & volTab$logFC > thresh, "diffCol"] <- "sigUp"
volTab[volTab$logFC < 0 & volTab$p.adj < 0.1 & volTab$logFC < -thresh, "diffCol"] <- "sigDown"
volTab[volTab$p.adj >= 0.1, "diffCol"] <- "ns"
volTab[volTab$logFC < thresh & volTab$logFC > -(thresh), "diffCol"] <- "ns"
p <- volTab %>%
mutate(logFC = ifelse(logFC > maxLim, maxLim, logFC)) %>%
mutate(logFC = ifelse(logFC < minLim, minLim, logFC)) %>%
#filter(gene %in% geneSelec) %>%
filter(Treatment == x) %>%
ggplot(aes(x = logFC, y = -log10(p.adj), fill = diffCol)) +
geom_point(shape = 21, size = 2) + xlab("Log2FC") + ylab("-Log10(p.adj)") + ggtitle(paste0(x)) +
#xlim(minLim, maxLim) +
geom_hline(yintercept = -log10(0.1), linetype = "dashed") +
geom_vline(xintercept = -thresh, linetype = "dashed") + geom_vline(xintercept = thresh, linetype = "dashed") +
geom_label_repel(data = filter(volTab, Symbol %in% c(geneSelec), Treatment %in% x, p.adj < 0.1,
logFC > thresh | logFC < -(thresh)), aes(label = Symbol), fill = "white",
nudge_y = 0.25, size = 4, max.overlaps = getOption("ggrepel.max.overlaps", default = 30),
segment.linetype = 2, force = 20, alpha = 0.8) +
scale_fill_manual(values = c("sigUp" = "#F44336", "sigDown" = "#1976D2", "ns" = "grey")) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
text = element_text(size = 17.5),
axis.text = element_text(size = 15),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
legend.background = element_rect(fill='transparent'),
strip.background = element_blank())
ggsave(plot = p, filename = paste0(opt$plot, "vol_", x, ".png"),
height = 15, width = 15, unit = "cm")
p
})
pList
## [[1]]
#### Plot bafilomycin A1 vs Motolimod
p <- de.res %>% filter(Treatment %in% c("BafilomycinA1", "Motolimod")) %>%
select(gene, Symbol, logFC, Treatment) %>%
pivot_wider(names_from = "Treatment", values_from = logFC) %>%
ggplot(aes(x = BafilomycinA1, y = Motolimod)) +
xlab("Bafilomycin A1 (log2FC)") + ylab("Motolimod (log2FC)") +
ggtitle("Motolimod vs Bafilomycin A1") +
geom_smooth(method = "lm", colour = "black", fill = "grey80") +
geom_point() +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
text = element_text(size = 17.5),
axis.text = element_text(size = 15),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
legend.background = element_rect(fill='transparent'),
strip.background = element_blank())
p
## `geom_smooth()` using formula = 'y ~ x'
ggsave(plot = p, filename = paste0(opt$plot, "moto_vs_bafilo.png"),
height = 15, width = 15, units = "cm")
## `geom_smooth()` using formula = 'y ~ x'
TP53 dot plot
geneSelec <- path.info %>% filter(pathway %in% c("TP53up")) %>% select(gene) %>% unlist() %>% as.vector()
geneSub <- lapply(unique(de.res$Treatment), function(x) {
de.res %>%
filter(Treatment == x) %>%
filter(p.adj < 0.1) %>%
filter(logFC > 0.25 | logFC < -0.25) %>%
pull(Symbol)
}) %>% unlist() %>% unique()
geneSub <- intersect(geneSelec, geneSub)
plotMat.wide <- de.res %>%
mutate(logFC = ifelse(p.adj >= 0.1, 0, logFC)) %>%
dplyr::select(Symbol, Treatment, logFC) %>%
filter(Symbol %in% geneSub) %>%
dplyr::rename(log2FC = logFC) %>%
pivot_wider(names_from = "Treatment", values_from = "log2FC") %>%
column_to_rownames("Symbol")
geneOrder <- rownames(plotMat.wide)[hclust(dist(plotMat.wide), method = "ward.D2")$order]
p <- de.res %>% filter(#Treatment %in% c("Birinapant", "Motolimod"),
Symbol %in% geneSub) %>%
#filter(Symbol == "MDM2") %>%
filter(p.adj < 0.1) %>%
mutate(pSign = -log10(p.adj) * sign(logFC),
p.adj = ifelse(-log10(p.adj) > 5, 5, -log10(p.adj)),
Symbol = factor(Symbol, levels = geneOrder)) %>%
dplyr::rename(log2FC = logFC) %>%
ggplot(aes(y = Symbol, x = Treatment)) +
geom_point(aes(size = p.adj, fill = log2FC), shape = 21) +
xlab("") + ylab("") + ggtitle("Effect on TP53 Target Genes") +
lgd +
scale_fill_gradient2(high = "#F44336", low = "#1976D2", name = "log2FC",
limits=c(-2, 2), oob=squish) +
scale_size_continuous(name = "-Log10(p.adj)",
breaks = seq(0, 5, by = 2.5),
limits = c(0, 5)
) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
text = element_text(size = 12.5),
panel.border = element_rect(linewidth = 1))
p
ggsave(plot = p, filename = paste0(opt$plot, "tp53_long.png"),
height = 40, width = 12.5, units = "cm")
GSEA
go.df.merge <- mclapply(c("Birinapant", "Motolimod", "Selinexor", "BafilomycinA1", "Dacinostat", "Pomalidomide", "Ibrutinib",
"Everolimus", "Venetoclax", "Ruxolitinib", "Nutlin-3a"), mc.cores = ncores,
function(compTreat) {
message(paste("Running gprofiler for", compTreat, sep = " "))
sigUp <- de.res %>% filter(logFC > 0.25, Treatment == compTreat,
#!gene %in% c(drugOther),
p.adj <= 0.1) %>% arrange(desc(logFC))
#sigUp <- de.res %>% filter(celltype == compType, logFC > 0.25, Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
sigUp
if(nrow(sigUp) > 1) {
gostresUp <- gost(query = sigUp$gene,
organism = "hsapiens", ordered_query = TRUE,
multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.1, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = NULL, as_short_link = FALSE)
gostresUp$result
go.df.up <- gostresUp$result %>% arrange(p_value)
go.df.up$rank <- 1:nrow(go.df.up)
go.df.up$type <- "Up"
} else {
go.df.up <- data.frame(p_value = NA, rank = NA, type = "Up", source = NA, term_name = NA,
intersection_size = 10)
}
sigDown <- de.res %>% filter(logFC < -(0.25), Treatment == compTreat,
#!gene %in% c(drugOther),
p.adj <= 0.1) %>% arrange(logFC)
#sigDown <- de.res %>% filter(celltype == compType, logFC < -(0.25), Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
if(nrow(sigDown) > 1) {
gostresDown <- gost(query = sigDown$gene,
organism = "hsapiens", ordered_query = TRUE,
multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.1, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = NULL, as_short_link = FALSE)
gostresDown$result
go.df.down <- gostresDown$result %>% arrange(p_value)
go.df.down$rank <- 1:nrow(go.df.down) * -1
go.df.down$type <- "Down"
} else {
go.df.down <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA,
intersection_size = 10)
}
go.df.merge <- rbind(go.df.up, go.df.down)
go.df.merge$Treatment <- compTreat
go.df.merge <- go.df.merge %>% filter(!intersection_size < 10) %>%
filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(p.adj = p.adjust(p_value, method = "BH"))
if(nrow(go.df.merge) < 1) {
go.df.merge <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA,
intersection_size = 0)
} else {
go.df.merge
}
}) %>% bind_rows()
## Exclude all terms that are up and downregulated by a single-drug
dblTerm <- go.df.merge %>%
filter(!intersection_size < 10) %>%
filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
#filter(term_name %in% unique(termTab$term_name)) %>%
dplyr::count(term_name, Treatment) %>%
filter(n > 1) %>%
pull(term_name) %>% unique()
plotMat <- go.df.merge %>%
mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
filter(!source %in% c("CORUM", "TF", "MIRNA"), !intersection_size < 10) %>%
filter(!term_name %in% c(dblTerm)) %>%
#filter(term_name %in% gseaSub) %>%
select(Treatment, term_name, type, p.adj) %>%
mutate(p.adj = -log10(p.adj),
p.adj = ifelse(type == "Down", p.adj * -1, p.adj)) %>%
select(-type) %>%
pivot_wider(names_from = "term_name", values_from = "p.adj") %>%
column_to_rownames("Treatment")
plotMat[is.na(plotMat)] <- 0
pHeat <- Heatmap(plotMat, show_column_names = FALSE,
clustering_method_rows = "ward.D2",
clustering_method_columns = "ward.D2", border = TRUE,
column_title = paste0("Enriched Gene Sets Among DEGs (FDR < 10%)"),
column_title_gp = gpar(fontsize = 17.5),
heatmap_legend_param = list(title_gp = gpar(fontsize = 12.5)),
col=colorRamp2(c(-5, 0, 5), colors = c("#1976D2", "white", "#F44336")),
name = "-Log10(p.adj)\nwith Direction")
## Warning: The input is a data frame-like object, convert it to a matrix.
pHeat
png(file=paste0(opt$plot, "GSEA_heatmap.png"),
height = 10, width = 25, res = 600, units = "cm", bg = "transparent")
draw(pHeat, background = "transparent")
dev.off()
## png
## 2
## Plot top enriched pathwaysgo.df.merge
nTop <- 20
lapply(c("Selinexor"), function(y) {
gseaTab <- go.df.merge %>% filter(Treatment == y) %>% filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
mutate(term_name = paste0(source, "_", term_name)) %>% unique()
message(y)
lapply(unique(gseaTab$type), function(z) {
message(z)
gseaTab <- gseaTab %>% filter(type == z)
gseaTab <- gseaTab[1:nTop, ] # get only top enrichments
termOrder <- gseaTab %>% arrange(desc(p_value)) %>% pull(term_name)
gseaTab <- gseaTab %>%
mutate(term_name = factor(term_name, levels = termOrder),
p.adj = p.adjust(p_value, method = "BH"))
if(-log10(min(gseaTab$p.adj, na.rm = TRUE)) > 10) {
x.max <- -log10(min(gseaTab$p.adj, na.rm = TRUE))
} else {
x.max <- 10
}
if(y %in% c("Birinapant", "Selinexor", "Pomalidomide")) {
x.max <- 3
}
if(y %in% c("Pomalidomide")) {
x.max <- 2
}
p <- gseaTab %>%
ggplot(aes(x = -log10(p.adj), y = term_name)) +
geom_bar(stat = "identity", colour = "black", alpha = 0.5, aes(fill = type)) +
geom_vline(xintercept = -log10(0.1), alpha = 0.8, linetype = "dashed") +
geom_text(aes(label = term_name, y = term_name, x = 0.25), hjust = 0, size = 4.25) +
xlab("-Log10(p.adj)") + ylab("") + #ggtitle(paste0("Top enriched gene-set: ", z)) +
ggtitle(paste0(y, " - ", z)) +
xlim(0, x.max) +
lgd + theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
text = element_text(size = 22.5),
panel.border = element_rect(colour = "black", fill=NA, size=1.5),
strip.background = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
legend.background = element_rect(fill='transparent'),
legend.position = "none") +
scale_fill_manual(values = c("Down" = "#1976D2", "Up" = "#F44336"))
#scale_fill_npg() +
#guides(fill=guide_legend(title="Source"))
ggsave(plot = p,
file = paste0(opt$plot, "GSEA/", y, "_", z, "_top_", nTop, ".png"),
height = 22.5, width = 19, units = "cm")
p
})
})
## [[1]]
## [[1]][[1]]
Hypergeometric test for transcription factors
treatSub <- c("Pomalidomide")
## Run hypergeometric test
hyper.df <- lapply(treatSub, function(x) {
## Limit the number of targets
motif.df <- motif.df %>% filter(n_targets < maxtargets)
message("Keeping ", nrow(motif.df), " TFs with less than ", maxtargets)
## Get the transcription factor targets
lapply(c("Up", "Down"), function(y) {
df.hub <- de.res %>% filter(p.adj < 0.05, Treatment == x)
if(y == "Up") {df.hub <- df.hub %>% filter(logFC > 0.25)}
if(y == "Down") {df.hub <- df.hub %>% filter(logFC < -0.25)}
mclapply(unique(motif.df$motif_name), mc.cores = ncores,
function(z) {
message(z)
tf.targets <- target.df %>% filter(motif_name == z) %>% pull(targets)
size.module <- nrow(df.hub)
size.genome <- length(unique(de.res$Symbol))
size.genome <- nrow(sce)
n.targets <- length(tf.targets)
n.over <- intersect(df.hub$Symbol, tf.targets) %>% length()
p.val <- phyper(q = n.over - 1, # number of white balls drawn
m = n.targets, # number white balls in urn
n = size.genome - n.targets, # number black balls in urn
k = size.module, # number of draws
lower.tail= FALSE)
df <- data.frame(Treatment = x, dir = y, tf = z, p.val = p.val, n.over = n.over,
n.targets = n.targets, size.module = size.module,
size.genome = size.genome)
#df %>% filter(n.over != 0)
}) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p.val, method = "BH"))
})# %>% bind_rows()
names(hyper.df) <- treatSub
lapply(c("Pomalidomide"), function(x) {
lapply(c("Up", "Down"), function(y) {
message(x)
hyper.df <- hyper.df %>% bind_rows() %>%
filter(Treatment == x, dir == y) %>%
arrange(p.adj) %>%
filter(p.adj < 0.1) %>%
mutate(p.adj = -log10(p.adj))
hyper.df$rank <- 1:nrow(hyper.df)
## Get the top 5 enriched TPs
toptf <- hyper.df %>% filter(rank <= 10) %>%
pull(tf)
if(x == "Pomalidomide") {
toptf <- c(toptf, "IKZF1", "IKFZ3", "BRD7", "ARID2", "IRF4", "MYC")
}
if(x %in% c("Nutlin-3a", "Selinexor")) {
toptf <- c(toptf, "TP53")
}
if(x %in% c("Birinapant", "Motolimod", "BafilomycinA1", "Ibrutinib")) {
toptf <- c(toptf, "NFKB", "NFKB1", "NFKB2", "REL", "RELA", "RELB")
}
if(y == "Up") {col.point <- "#F44336"}
if(y == "Down") {col.point <- "#1976D2"}
p <- hyper.df %>%
ggplot(aes(x = rank, y = p.adj)) +
geom_point() +
geom_point(data = hyper.df[hyper.df$tf %in% toptf, ], fill = col.point, size = 2.5, shape = 21) +
xlab("Rank") + ylab("-Log10(p.adj)") + ggtitle(paste0(x, " - TF Enrichment (", y, ")")) +
geom_label_repel(aes(label = tf), data = hyper.df[hyper.df$tf %in% toptf, ],
max.overlaps = getOption("ggrepel.max.overlaps", default = 30),
segment.linetype = 2, force = 20, alpha = 0.8) +
lgd +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
text = element_text(size = 17.5),
axis.text = element_text(size = 15),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
legend.background = element_rect(fill='transparent'),
strip.background = element_blank())
ggsave(plot = p, filename = paste0(opt$plot, x, "_", y, "_tf_enrichment.png"),
height = 12.5, width = 16, units = "cm")
p
})
})
## [[1]]
## [[1]][[1]]
##
## [[1]][[2]]
Overview scRNA-Seq data
## Randomly sample
set.seed(100)
sce.sc <- sce.sc[, sample(1:ncol(sce.sc), ncol(sce.sc))]
## Re-run UMAP
set.seed(100)
sce.sc <- runUMAP(sce.sc, n_neighbors = 15, min_dist = 0.75, name = "UMAP", BPPARAM = mcparam)
## Make plots
plotTab <- ggcells(sce.sc)[[1]]
## Pairwise each treatment
lapply(c("Birinapant", "Nutlin-3a", "Ibrutinib", "Everolimus", "Selumetinib"), function(x) {
p <- plotTab %>%
filter(Treatment %in% c("DMSO", x)) %>%
ggplot(aes(x = UMAP.1, y = UMAP.2, fill = Treatment)) +
geom_point(shape = 21, size = 2) +
lgd + theme_void() + theme(legend.position = "none") +
scale_fill_manual(values = c("DMSO" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000",
"Ibrutinib" = "#283593", "Everolimus" = "#43A047",
"Selumetinib" = "#F8BBD0"))
ggsave(plot = p, filename = paste0(opt$plot, x, "_scrna_umap.png"),
height = 15, width = 15, units = "cm")
p
})
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Show donut plots
plotTab.vdj <- plotTab %>%
left_join(types, by = "uniqueBarcode") %>%
filter(celltype %in% c("malignant B-cell", "healthy B-cell")) %>%
mutate(celltype = ifelse(celltype == "malignant B-cell", patientID, celltype),
celltype = ifelse(celltype == "healthy B-cell", "Healthy B-Cells", celltype))
plotTab.vdj[plotTab.vdj$celltype == "H525", "celltype"] <- "HCL1"
plotTab.vdj[plotTab.vdj$celltype == "H526", "celltype"] <- "HCL2"
plotTab.vdj[plotTab.vdj$celltype == "H432", "celltype"] <- "MCL1"
plotTab.vdj[plotTab.vdj$celltype == "H452", "celltype"] <- "MCL2"
plotTab.vdj[plotTab.vdj$celltype == "H496", "celltype"] <- "MCL3"
## Merge VDJ data with plotting info
plotTab.vdj <- plotTab.vdj %>% left_join(select(immTab, cell_id, clone_id, c_call#, mu_freq
),
by = c("uniqueBarcode" = "cell_id"))
plotTab.vdj %>% filter(clone_id %in% c("60_5")) %>% select(celltype)
## [1] celltype
## <0 rows> (or 0-length row.names)
## Calculate the number of clones per cell type
countTab <- plotTab.vdj %>%
group_by(celltype) %>%
select(celltype, clone_id) %>% unique() %>%
dplyr::count(celltype) %>% mutate(n_clones = n) %>%
select(-n)
countTab
## # A tibble: 6 × 2
## # Groups: celltype [6]
## celltype n_clones
## <chr> <int>
## 1 HCL1 2
## 2 HCL2 2
## 3 Healthy B-Cells 176
## 4 MCL1 2
## 5 MCL2 2
## 6 MCL3 2
## Calculate the size of each clone
countSize <- plotTab.vdj %>%
filter(!is.na(clone_id)) %>%
group_by(celltype) %>%
select(celltype, clone_id) %>%
dplyr::count(celltype, clone_id) %>% mutate(clone_size = n) %>%
select(-n)
countSize
## # A tibble: 180 × 3
## # Groups: celltype [6]
## celltype clone_id clone_size
## <chr> <chr> <int>
## 1 HCL1 BCR_122_142 3006
## 2 HCL2 BCR_181_7 2564
## 3 Healthy B-Cells BCR_100_64 1
## 4 Healthy B-Cells BCR_101_82 1
## 5 Healthy B-Cells BCR_102_135 1
## 6 Healthy B-Cells BCR_103_129 1
## 7 Healthy B-Cells BCR_104_54 1
## 8 Healthy B-Cells BCR_105_46 1
## 9 Healthy B-Cells BCR_106_163 1
## 10 Healthy B-Cells BCR_107_34 1
## # ℹ 170 more rows
###################
# Donut charts of clonal composition
###################
## Calculate the size of each clone - don't subset to only cells with clone_id
# otherwise NA T-cells will not be in the plot later.
cloneCount <- plotTab.vdj %>%
group_by(celltype) %>%
dplyr::count(clone_id) %>% ungroup()
cloneCount
## # A tibble: 186 × 3
## celltype clone_id n
## <chr> <chr> <int>
## 1 HCL1 BCR_122_142 3006
## 2 HCL1 <NA> 626
## 3 HCL2 BCR_181_7 2564
## 4 HCL2 <NA> 513
## 5 Healthy B-Cells BCR_100_64 1
## 6 Healthy B-Cells BCR_101_82 1
## 7 Healthy B-Cells BCR_102_135 1
## 8 Healthy B-Cells BCR_103_129 1
## 9 Healthy B-Cells BCR_104_54 1
## 10 Healthy B-Cells BCR_105_46 1
## # ℹ 176 more rows
cloneNum <- lapply(unique(plotTab.vdj$celltype), function(x) {
subTab <- filter(plotTab.vdj, celltype == x) # subset to each patient
countTab <- dplyr::count(subTab, clone_id)
nclones <- countTab %>% nrow() # get the number of clones
ncells_bcr <- nrow(subTab) # calculate the number of B or T cells per patient - important to subset plotTab to only B or T-cell clusters before doing this
res <- data.frame(patID = x, nclones = nclones, ncells_bcr = ncells_bcr)
res$type <- opt$subType
res
}) %>% bind_rows()
cloneNum
## patID nclones ncells_bcr
## 1 MCL1 2 3084
## 2 HCL1 2 3632
## 3 HCL2 2 3077
## 4 MCL3 2 3799
## 5 MCL2 2 2082
## 6 Healthy B-Cells 176 218
df <- left_join(cloneCount, cloneNum, by = c("celltype" = "patID")) %>%
group_by(celltype) %>% mutate(clone_ratio = n / ncells_bcr) %>% # divide the size of the clone
select(celltype, clone_id, clone_ratio, nclones) %>% unique() %>%
mutate(ymax = cumsum(clone_ratio), ymin = c(0, head(ymax, n= -1)))
df
## # A tibble: 186 × 6
## # Groups: celltype [6]
## celltype clone_id clone_ratio nclones ymax ymin
## <chr> <chr> <dbl> <int> <dbl> <dbl>
## 1 HCL1 BCR_122_142 0.828 2 0.828 0
## 2 HCL1 <NA> 0.172 2 1 0.828
## 3 HCL2 BCR_181_7 0.833 2 0.833 0
## 4 HCL2 <NA> 0.167 2 1 0.833
## 5 Healthy B-Cells BCR_100_64 0.00459 176 0.00459 0
## 6 Healthy B-Cells BCR_101_82 0.00459 176 0.00917 0.00459
## 7 Healthy B-Cells BCR_102_135 0.00459 176 0.0138 0.00917
## 8 Healthy B-Cells BCR_103_129 0.00459 176 0.0183 0.0138
## 9 Healthy B-Cells BCR_104_54 0.00459 176 0.0229 0.0183
## 10 Healthy B-Cells BCR_105_46 0.00459 176 0.0275 0.0229
## # ℹ 176 more rows
#### Overview over all patients/celltypes
pClone.Pat <- df %>%
ggplot(aes(ymax = ymax, ymin = ymin, xmax=4, xmin=3, fill = clone_id)) +
geom_rect(colour = "black", size = 0.25) +
coord_polar(theta="y") +
xlim(c(2, 4)) +
# ggtitle(paste0("Clonal Composition per Patient - ", opt$subType)) +
facet_wrap(. ~ celltype, nrow = 1) +
theme_void() +
theme(plot.title = element_text(hjust = 0.5),
strip.background = element_blank(),
strip.text.x = element_text(size = 12.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
legend.position = "none")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
pClone.Pat
#### Output individual patients
subPat <- c("HCL1", "HCL2", "MCL1", "MCL2", "MCL3", "Healthy B-Cells")
lapply(subPat, function(subPat) {
pClone <- df %>% filter(celltype == subPat) %>%
ggplot(aes(ymax = ymax, ymin = ymin, xmax=4, xmin=3, fill = clone_id)) +
geom_rect(colour = "black", size = 0.25) +
coord_polar(theta="y") +
xlim(c(2, 4)) +
ggtitle(paste0(subPat)) +
theme_void() +
theme(plot.title = element_text(hjust = 0.5, size = 75),
legend.position = "none")
pClone
ggsave(plot = pClone, filename = paste0(opt$plot, "bcr_comp_", subPat, ".png"),
height = 20, width = 20, units = "cm")
})
## [[1]]
## [1] "plots/SFig4_5/bcr_comp_HCL1.png"
##
## [[2]]
## [1] "plots/SFig4_5/bcr_comp_HCL2.png"
##
## [[3]]
## [1] "plots/SFig4_5/bcr_comp_MCL1.png"
##
## [[4]]
## [1] "plots/SFig4_5/bcr_comp_MCL2.png"
##
## [[5]]
## [1] "plots/SFig4_5/bcr_comp_MCL3.png"
##
## [[6]]
## [1] "plots/SFig4_5/bcr_comp_Healthy B-Cells.png"
Per patient UMAPs
patSub <- names(sce.pat)
sce.pat <- mclapply(patSub, mc.cores = ncores, function(x) {
sce <- sce.pat[[x]]
set.seed(100)
sce <- runUMAP(sce, n_neighbors = 15, min_dist = 0.35, name = "UMAP", BPPARAM = mcparam)
})
names(sce.pat) <- patSub
patSub <- c("H371", "H431", "H279")
lapply(patSub, function(x) {
sce <- sce.pat[[x]]
set.seed(100)
sce <- sce[,sample(1:ncol(sce), ncol(sce))]
plotTab <- ggcells(sce)[[1]]
if(x == "H371") {tlt <- "T-PLL1"
} else if(x == "H431") {tlt <- "T-PLL2"
} else if(x == "H279") {tlt <- "T-PLL3"} else {
tlt <- x
}
p <- plotTab %>%
ggplot(aes(x = UMAP.1, y = UMAP.2, fill = Treatment)) +
geom_point(shape = 21, size = 3) +
ggtitle(tlt) +
lgd +
theme_void() +
theme(axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = "none",
plot.title = element_text(size = 20, hjust = 0.5)) +
scale_fill_manual(values = c("DMSO" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000",
"Ibrutinib" = "#283593", "Everolimus" = "#43A047",
"Selumetinib" = "#F8BBD0"))
ggsave(plot = p, filename = paste0(opt$plot, x, "_perpat_umap.png"),
height = 12.5, width = 12.5, units = "cm")
p
})
## [[1]]
##
## [[2]]
##
## [[3]]
Output session info
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel grid stats4 stats graphics grDevices datasets
## [8] utils methods base
##
## other attached packages:
## [1] sna_2.7-2 statnet.common_4.9.0
## [3] scales_1.3.0 numbat_1.4.0
## [5] Matrix_1.6-1.1 airr_1.5.0
## [7] ggsci_3.2.0 BiocParallel_1.36.0
## [9] ggnet_0.1.0 network_1.18.2
## [11] tftargets_1.3 ggrepel_0.9.5
## [13] readxl_1.4.3 edgeR_4.0.16
## [15] limma_3.58.1 circlize_0.4.16
## [17] ComplexHeatmap_2.18.0 gridExtra_2.3
## [19] gprofiler2_0.2.3 scater_1.30.1
## [21] scran_1.30.2 scuttle_1.12.0
## [23] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
## [25] Biobase_2.62.0 GenomicRanges_1.54.1
## [27] GenomeInfoDb_1.38.8 IRanges_2.36.0
## [29] S4Vectors_0.40.2 BiocGenerics_0.48.1
## [31] MatrixGenerics_1.14.0 matrixStats_1.3.0
## [33] lubridate_1.9.3 forcats_1.0.0
## [35] stringr_1.5.1 dplyr_1.1.4
## [37] purrr_1.0.2 readr_2.1.5
## [39] tidyr_1.3.1 tibble_3.2.1
## [41] ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.22 splines_4.3.2
## [3] bitops_1.0-7 ggplotify_0.1.2
## [5] cellranger_1.1.0 polyclip_1.10-7
## [7] lifecycle_1.0.4 doParallel_1.0.17
## [9] lattice_0.21-9 MASS_7.3-60
## [11] magrittr_2.0.3 plotly_4.10.4
## [13] sass_0.4.9 rmarkdown_2.27
## [15] jquerylib_0.1.4 yaml_2.3.9
## [17] metapod_1.10.1 RColorBrewer_1.1-3
## [19] abind_1.4-5 zlibbioc_1.48.2
## [21] quadprog_1.5-8 hahmmr_1.0.0
## [23] ggraph_2.2.1 RCurl_1.98-1.14
## [25] yulab.utils_0.1.5 tweenr_2.0.3
## [27] GenomeInfoDbData_1.2.11 irlba_2.3.5.1
## [29] tidytree_0.4.6 dqrng_0.4.1
## [31] DelayedMatrixStats_1.24.0 codetools_0.2-19
## [33] DelayedArray_0.28.0 ggforce_0.4.2
## [35] tidyselect_1.2.1 shape_1.4.6.1
## [37] aplot_0.2.3 farver_2.1.2
## [39] ScaledMatrix_1.10.0 viridis_0.6.5
## [41] jsonlite_1.8.8 GetoptLong_1.0.5
## [43] BiocNeighbors_1.20.2 tidygraph_1.3.1
## [45] iterators_1.0.14 systemfonts_1.1.0
## [47] foreach_1.5.2 tools_4.3.2
## [49] ragg_1.3.2 treeio_1.29.0.002
## [51] Rcpp_1.0.12 glue_1.7.0
## [53] SparseArray_1.2.4 mgcv_1.9-0
## [55] xfun_0.45 withr_3.0.0
## [57] fastmap_1.2.0 bluster_1.12.0
## [59] fansi_1.0.6 digest_0.6.36
## [61] rsvd_1.0.5 parallelDist_0.2.6
## [63] timechange_0.3.0 R6_2.5.1
## [65] gridGraphics_0.5-1 textshaping_0.4.0
## [67] colorspace_2.1-0 Cairo_1.6-2
## [69] RhpcBLASctl_0.23-42 utf8_1.2.4
## [71] generics_0.1.3 renv_1.0.7
## [73] data.table_1.15.4 graphlayouts_1.1.1
## [75] httr_1.4.7 htmlwidgets_1.6.4
## [77] S4Arrays_1.2.1 uwot_0.2.2
## [79] pkgconfig_2.0.3 gtable_0.3.5
## [81] XVector_0.42.0 htmltools_0.5.8.1
## [83] clue_0.3-65 png_0.1-8
## [85] ggfun_0.1.5 knitr_1.48
## [87] rstudioapi_0.16.0 tzdb_0.4.0
## [89] rjson_0.2.21 coda_0.19-4.1
## [91] nlme_3.1-163 cachem_1.1.0
## [93] GlobalOptions_0.1.2 vipor_0.4.7
## [95] pillar_1.9.0 logger_0.3.0
## [97] vctrs_0.6.5 BiocSingular_1.18.0
## [99] beachmat_2.18.1 cluster_2.1.4
## [101] beeswarm_0.4.0 evaluate_0.24.0
## [103] cli_3.6.3 locfit_1.5-9.10
## [105] compiler_4.3.2 rlang_1.1.4
## [107] crayon_1.5.3 labeling_0.4.3
## [109] fs_1.6.4 ggbeeswarm_0.7.2
## [111] stringi_1.8.4 viridisLite_0.4.2
## [113] munsell_0.5.1 lazyeval_0.2.2
## [115] hms_1.1.3 patchwork_1.2.0
## [117] sparseMatrixStats_1.14.0 statmod_1.5.0
## [119] highr_0.11 igraph_2.0.3
## [121] memoise_2.0.1 scistreer_1.2.0
## [123] RcppParallel_5.1.8 bslib_0.7.0
## [125] phangorn_2.11.1 ggtree_3.10.1
## [127] fastmatch_1.1-4 ape_5.8